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ABSTRACT 

The multicanonical ansatz is used to study variations in the energy landscape of a small 
peptide, Met-enkephalin, under a change from the ECEPP/2 force field to ECEPP/3. Lo- 
cal minima with energies up to 5 kcal/mol higher than the global minima are sampled 
and classified according to H-bridges and backbone angles. The distribution and relative 
weight for various temperatures of the minima are calculated and compared for the two 
force fields. We demonstrate that while there are small differences in the energy land- 
scape our results at relevant temperatures are robust under changes between ECEPP/2 
to ECEPP/3. 
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Introduction 



While there is considerable progress in numerical simulation of peptides and proteins (for 
a recent review see, for instance, 0), prediction of their low temperature conformations 
solely from first principles remains a formidable task. First, the interactions between 
the atoms in the molecules are described by force fields which depend in part on empir- 
ically determined parameters. It is still an open question which of the many force fields 
(AMBER §, CHARMM §, ECEPP § - [0, ...) is the optimal choice (or if they are 
sufficiently accurate at all). The other mayor problem is common to many important 
physical systems. At low temperatures the energy landscape of proteins is characterized 
by a multitude of local minima separated by high energy barriers. In simulations based on 
classical molecular dynamics or Monte Carlo techniques || these barriers can be seldom 
crossed. Hence, only small parts of phase space are sampled (in a finite number of Monte 
Carlo sweeps) and physical quantities cannot be calculated accurately. For this reason, 
the "native" conformation of a protein (which should correspond to the global minimum 
in free energy) is often identified with the lowest potential energy configuration. In this 
approximation entropic contributions are ignored, but it allows application of optimiza- 



tion techniques like simulated annealing ||10|| , genetic algorithms |9[ or Monte Carlo with 
Minimization to tackle the protein folding problem. However, this approximation 
may be to crude. Especially for smaller peptides entropic contributions are expected to 
be important. 

It is not always clear whether the limited quality of numerical results is due to insuffi- 
cient simulation algorithms or inadequate force fields. Recent progress in the development 
of Monte Carlo techniques may change this situation. Simulations in generalized ensem- 
bles promise a much better sampling of the phase space. A numerical comparison of some 
of these new algorithms can be found in Ref. |12|]. One of its more prominent expo- 
nents is the so called multicanonical approach of Berg and co-workers |R| who propose 
a weighting scheme which yield to a flat probability distribution in energy. Hence, all 
energies have equal weight and an one-dimensional random walk in energy space is car- 
ried out (when simulated with local updates) which ensures that the simulation will not 
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get trapped in any local energy minimum. The Boltzmann distributions may be obtained 
for a given range of temperatures from one multicanonical simulation by re-weighting the 



states [14]. The method was exploited to model first-order phase transitions [13, 15] and 



spin-glass systems fTR O, |Ig|. The prediction of peptide and protein three-dimensional 



structures with multicanonical algorithms was first addressed in Ref. [19| for Monte Carlo 



methods and in Ref. [20] for molecular dynamics. Subsequent works include the study 



of coil-globular transitions of a model protein [|2l|, helix-coil transitions of amino acid 
homo-oligomers |22| and the conformational sampling of a constrained peptide [^3j . 



In the current article we use the multicanonical technique to compare distribution and 
relative weight of the local minima at room temperature and below for the closely related 
ECEPP/2 and ECEPP/3 force field in the case of the small pentapeptide Met-enkephalin. 
ECEPP/3 deviates from the previous ECEPP/2 by slightly different parameters due to 
improvements in the experimental measurements. We want to study the influence of 
such small changes on the energy landscape and our observed quantities. This should 
show us how much numerical simulations are affected by experimental uncertainties in 
the parameters of the force fields. 

Methods 

Multicanonical Algorithm 

Although the multicanonical approach is explained in detail elsewhere (see, for instance 
Refs. |T6| , p4]]), we briefly summarize the idea and implementation of the method for 
completeness. 

Simulations in the canonical ensemble weight each configuration with wb{E, j3) = e~" E 
and yield the usual Boltzmann probability density distribution of energies: 

P B (E,T)ocn(E)e-? E (1) 

where n(E) is the density of states with energy E, and (3 = with temperature 
T and Boltzmann constant ks- 

On the other hand, in the multicanonical approach configurations with energy E are 
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updated with a weight: 

w mu (E) <xn-\E) = e' s ^ (2) 
where S(E) is the microcanonical entropy. A uniform distribution of the energy 

P mu {E) = n(E)w{E) = const (3) 

is obtained from the simulation with these weight factors. In this way information is 
collected over the whole energy range and from a single simulation one can obtain the 
canonical distribution for a wide range of temperatures by the re- weighting techniques: |14 



P B (T,E)cxP mu (E)w-ie^ E . (4) 



This allows one to calculate any thermodynamic quantity O by 



dE 0(E)P B (T,E) 

<0> T = J — r . (5) 

J dE P B (T,E) 

However, unlike in a canonical simulation the weights w mu (E) are not a priori known 
(in fact, knowledge of the exact weights is equivalent with solving the system), and one 
needs its estimators for a numerical simulation. One way to obtain these estimators is 
the following iterative procedure. Starting with an initial guess w^^E) = e~ l3 ° E ', i.e., 
performing a canonical simulation at sufficient high temperature T , improved estima- 
tors of the multicanonical weights are calculated from histograms P'^}{E) of preceding 
simulations by 

w mu {E) = pi _^ E y (6) 



For details, see Ref. P21 and [124] . This method for calculating multicanonical weights is 



by no means unique and while it is quite general, it has the disadvantage that it requires 
a certain number of extra iterations which is not a priori known. 

Force-Fields 

In ECEPP the potential energy function E tot is given by the sum of the electrostatic 

term Ec, 12-6 Lennard- Jones term Elj, and hydrogen-bond term E BB for all pairs of 
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atoms in the peptide together with the torsion term E tor for all torsion angles: 

E to t = E c + E LJ + E HB + E tor , (7) 
Ec = (8) 

^ = (9) 

^ = eS-^). (io) 

(ij) vy '*i / 

E tor = J2 U i( 1±C0 < n iXi))- (11) 
i 

Here, is the distance between the atoms i and j, and xz is the torsion angle for the 
chemical bond I. In ECEPP bond lengths and bond angles (which are hard degrees of free- 
dom) are fixed at experimental values and no out-of-plane deformation of peptide bonds 
is allowed leaving the dihedral angles <f), if), u and x as independent variables. The various 
parameters A^, Bij, C^, Dij, C/j, and n{) for the energy function were determined by 
a combination of a priori calculations and minimization of the potential energies of the 
crystal lattices of single amino acids. In the light of more recent experimental findings the 
standard geometry and some energy parameters for prolyl and hydroxyprolyl have been 
updated. Together with a re-calculation of partial atomic charges of mainchain atoms 
this led to a revision of the ECEPP/2 parameter set to ECEPP/3. In the original formu- 
lation of the ECEPP/2 parameter set, the net charges of both the N- and the C-terminal 
residues, added up to the desired net charge of the molecule. Now each terminal amino 
acid residue carries a total charge of 0/+1/-1 at the N- or C-terminus, respectively. The 
correction avoids possible artifacts for systems with terminal groups of different types 
(e.g. one terminus charged, the other uncharged). 

Simulation and Technical Details 

To investigate changes in the energy landscape under small variations of the force field 
we have studied one of the simplest peptides, Met-enkephalin which has the amino acid 
sequence Tyr-Gly-Gly-Phe-Met. This peptide is convenient for our purpose, since for the 



potential energy function ECEPP/2 the lowest energy conformation is known |TT|, |25], ^| 
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and local minima with energies not much higher than the global minimum were sampled 



and classified for this molecule by Braun and co-workers. ||27|| We compare our results 
with theirs and calculate in addition the relative weight of the local minima at various 
temperatures which is not possible with the method used by Braun et al. 

For our simulations the backbone was terminated by a neutral -NH2 group at the 
N-terminus and a neutral -COOH group at the C-terminus as in the previous works of 
Met-enkephalin. For Met-enkephalin the two versions of ECEPP differ in the point charges 
for terminal groups which are listed in Tab. 1. The peptide-bond dihedral angles u were 
fixed at the value 180°, leaving 19 dihedral angles as independent variables.0 Interactions 
of the molecule with solvent molecules were neglected and the dielectric constant e was 
set equal to 2. We used the program SMC |28| which was modified to accommodate the 
multicanonical ensemble. 

The multicanonical weight factors were determined separately for both ECEPP/2 and 
ECEPP/3 by the iterative procedure described above. We needed 100,000 sweeps for 
their calculation. One MC sweep updates every dihedral angle of the molecule once. 
For each update of an angle a new randomly chosen value out of the interval [— n, 71] 
was proposed and the new configuration then accepted or rejected by the Metropolis 
criteria. All thermodynamic quantities were calculated for ECEPP/2 and ECEPP/3 
separately from a production run of 200,000 MC sweeps following additional 10,000 sweeps 
for equilibration each. The energy after each sweep was stored for further analysis. In 
all cases, each simulation started from a completely random initial conformation ("Hot 
Start"), but in the case of ECEPP/2 we also checked our results for shorter runs against 
these with initial configuration of the known groundstate ("Cold Start"), and we found 
that the results are in agreement with those from random initial conformations. This 
suggests that thermal equilibrium has been attained in our simulations. This conjecture 
is further supported by Fig. 1, in which we show the timeseries of the energy for both 



2 We fixed these hard degrees of freedom to reduce the number of variables. It would be more rigorous 
to allow for bending of the peptide angles. Similarly, one could argue that the fixed geometry of ECEPP 
should be replaced by a flexible one (as in other force fields), for instance by including Fixman potentials. 
However, in this article we are mostly interested in comparing two variants of a force field, not in 
reproducing experimentally found structures by numerical simulation, which allow us to choose a less 
rigorous approach. 
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production runs. In the course of the simulation, a Id random walk in the energy between 
low energy states and high energy states is performed as one would expect for local updates 
from the definition of the multicanonical ensemble. It is evident that low energy states 
which are separated in the timeseries by high energy states (which correspond to high 
temperature states) are uncorrelated. The number of such "tunnel" events is therefore a 
lower limit for the number of independent low energy states visited in the simulation. 

In order to investigate the energy landscape of our peptide we not only stored the 
energy of the actual conformation after each MC sweep, but also minimized a copy of 
this conformation if its energy was less than a certain limit, kcal/mol. While this limit 
is somehow arbitrary, our experience shows that it allows a rough distinction between a 
"low temperature region" and a "high temperature region" (it is the expectation value of 
Energy at T = 300/r where the specific heat has its maximum). From the minimization 
process we receive a local minimum conformation. Configurations which yield the same 
local minima belong to the same valley in the energy landscape and are therefore related 
to each other. Since we are mainly interested in local minima not too far away from the 
global minimum and to save disc space, we introduced a second limit. Only if the energy 
of the minimized configuration was below —6 kcal/mol (~ 5 kcal/mol above the global 
minimum) it was stored for further analysis. In Tab. 2 we summarize the absolute number 
of local minima sampled in this way for both variants of the ECEPP force field. 

While the set of configurations we received from our production runs are all local 
minima, due to limitations of the minimizer (a Newton-Raphson variant) they are not 
necessarily the lowest configurations associated with a certain valley in the energy land- 
scape. Each valley is in itself rough and consists of subvalleys separated by energy barriers. 
Since we are not interested in the microscopic details of the energy landscape we further 
clustered our local minima and identified each cluster by its member with the lowest en- 
ergy. In this way, only a small number (less than 5) of groups remains at the end of our 
clustering by which the overwhelming number of local minima is classified. 

We describe now briefly our clustering technique. In a first step we sorted all minimized 
conformers according to their energy. In a second step, we looked for identical conformers. 
Two conformers were regarded as identical if they differ in no angle more than 1° (taking 
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symmetries in account) and in their energies less than 0.0001 kcal/mol. For further 
clustering, we kept only the non-identical conformations. Their number which is also 
shown in Tab. 2 is our estimate for the total number of local minima with energy below 
—6 kcal/mol. It is interesting to observe that the number of non-identical conformers 
sampled in the multicanonical runs is similar for each choice of force field while the 
absolute number of sampled local minima varies much more (see Tab. 2). In a third 
step, using all backbone angles, the conformer (from our remaining set of non-identical 
local minima conformations) with the lowest energy is compared to all other with higher 
energies, according to differences in their backbone dihedrals. The conformers having all 
their angles within 10° in comparison with the first one, with respect to symmetry, are 
put into the same cluster, which is represented by the first conformer. This procedure 
is repeated for the next conformer with higher energy, which was not counted in the 
first cluster by comparing its angles with all the other remaining "free" conformers. In 
this way, we are left with 160 clusters for ECEPP/2 and 571 clusters for ECEPP/3. In 
the final step these clusters were classified into few groups using similarities in hydrogen 
bonding and the significant dihedrals, necessary to preserve the characteristic hydrogen 
bonds. The final groups which should represent the main valleys in the energy landscape 
of our peptide are characterized by their lowest-energy conformer. 

In Fig. 2a and 2b we display how in the course of our simulations different valleys of 
the energy landscape (characterized by its corresponding group number) are visited. We 
show only that times when the simulation visit one of the 4 (3) most important clusters 
for ECEPP/2 (ECEPP/3) and no plot symbol means that the simulation did not visit 
one of these valleys. As one can see there are different slices in which the multicanonical 
simulation moved only inside of a valley separated by jumps between different valleys. 
Again the number of these slices will give a lower limit on the number of independent low 
energy states visited in the simulation. 

Results and Discussion 

Using the techniques described above we were able for each choice of the force field to 
classify the large number of local minimum conformations sampled in our multicanonical 
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runs into few groups which represent more than 80% of the sampled local minimum 
conformations and which correspond to the important valleys in the energy landscape of 
Met-enkephalin. In Fig. 3a and Fig. 3b the dominant groups for both variants of the 
ECEPP force field are shown. In addition, we present in Tab. 3 the dihedral angles for 
the global minimum conformations obtained by our method. These structures are also 
shown in Fig. 4a and Fig. 4b. 

The global minimum for ECEPP/3 has characteristic hydrogen bond of the NH-group 
of Phe(4) with the carbonyl group of Try(l) resulting in a f3 -bend between the first 
and fourth residue. To preserve this hydrogen bond only the backbone angles of Gly(2) 
and Gly(3) have to be conserved and therefore one finds a large number of variants 
differing from the global minimum stucture in the other backbone and sidechain dihedrals 
in this group (which we call group "1"). The second group ("2") is characterized by the 
local minimum configuration with the second lowest energy in the ECEPP/3 force field. 
This configuration has two possible hydrogen bonds between N-H and C=0 groups of 
the Gly(2) and Met (5) residue forming a H'-type /3-turn between the second and last 
residue. Since more backbone dihedrals have to be conserved to preserve the two hydrogen 
bonds the configurations of this group vary much less than in the case of the previous 
group. Despite their differences (see Figs. 3 and 4) the lowest energy exponent of both 
groups differs by only 0.04 kcal/mol in potential energy, so that one can regard the 
global minimum of Met-enkephalin in the ECEPP/3 force field as degenerated. However, 
the electrostatic energy for the ECEPP/3-minimum is raised substantially if ECEPP/2 
charges are used for this structure. The total potential energy increases from E tot = 
-10.84 kcal/mol (ECEPP/3) to E tot = -9.46 kcal/mol (ECEPP2) while the energy 
of the other configuration changes only from E tot = —10.81 kcal/mol (ECEPP/3) to 
Etot = —10.72 kcal/mol (ECEPP/2), which is the global minimum of Met-enkephalin in 
the ECEPP/2 force field |2§ Configurations of both group "1" and "2" were also 
found in our ECEPP/2 simulation. 

We remark that our peptide seems to have a much richer spectrum of low energy 
local minimum configurations in the ECEPP/3 force field than in the case of ECEPP/2, 
but show in Fig. 3b only one more (the other contribute with less than 1% each) group. 
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Conformers of this group ("H") are characterized by a hydrogen bond between C=0 
group of the Tyr(l) and the N-H group of Met (5) residues. Together with the backbone 
dihedrals of Gly(2), Gly(3) and Phe(4) this hydrogen bonding resembles that of an a- 
helix. Its lowest energy exponent has a potential energy of —10.56 kcal/mol which is 
only 0.25 kcal/mol higher than the groundstate. No such conformation was found in 
the ECEPP/2 simulation. In the case of ECEPP/2, two more groups appeared with a 
frequency of more than 1% and are shown in Fig. 3a. Both are IF type /5-turns. Group 
"a" is characterized by single hydrogen bond between Gly(2) and Met (5). The lowest 
energy conformation in this group has a potential energy of —9.74 kcal/mol which make 
this group to the one with the second lowest potential energy, separated by ~ 1 kcal/mol 
from the groundstate. Group "b" has its hydrogen bond between Tyr(l) and Phe(4), its 
lowest energy exponent has a potential energy of —9.52 kcal/mol). Both groups were also 
found in the ECEPP/3 simulation, but with frequencies less than 1% of the non-identical 
conformers. 

For ECEPP/2 our local minimum conformations resemble those found in previous 
work |27| where in difference to the present work the peptide angles were released. Fixing 
these angles therefore seems to be a reasonable approximation. But while our classification 
of low energy local minima can also be obtained with other methods, the multicanonical 
approach allows in addition the estimation of the relative weight of these conformers as 
a function of temperature, i.e. entropic contributions can be considered. As an example 
we show in Fig. 5a and 5b the relative weight of the dominant groups of conformation for 
both ECEPP/2 and ECEPP/3. Note that around 300-fT the percentage of configurations 
which are in none of our groups (and hence do not belong to the important valleys) 
decreases dramatically for both force fields. At high temperatures coil structures dominate 
which show large flexibility and and therefore a broad energy spectrum. On the other 
hand, at low temperatures ordered structures are expected which are each confined to 
one of the valleys in the energy landscape. In the case of the ECEPP/2 simulation 
the percentage of configurations of type '1' (the ECEPP/3 groundstate) never exceeded 
4%. Configurations of type '2' (the groundstate for ECEPP/2) are found with much 
higher probability, but at room temperature they contribute to only about 30 % of the 
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conformers which confirms older work [jT9|j . For ECEPP/3 we find that the probability 
of configurations of type '2' is comparable to that of the ECEPP/2 simulation. At room 
temperature they contribute again to around 30% of the conformers. Conformations of 
type '1' appear with much higher probability than in the ECEPP/2 simulation, but that of 
type '2' are still dominating in the range of shown temperatures. At relevant temperatures 
the probability for finding the ECEPP/2 groundstate (type '2') is twice as high as that 
of finding the ECEPP/3 groundstate (type '1'). Of course we expect that at T — K 
we will find only the ECEPP/3 groundstate (type '1'), however the energy difference 
is much too small (about 0.04 kcal/mol) to favor them at finite temperatures against 
type '2' structures. This proofs that at least for small peptides entropic contributions 
cannot be neglected. Any pure optimization method (in the potential energy) may lead 
to configurations which are of only limited significance at relevant temperatures. 

In the final plots we like to show how the variations in the force field affects two 
important thermodynamic quantities. Fig. 6 displays the average energy < E > as a 
function of temperature. Again, we observe only little differences between the ECEPP/2 
and ECEPP/3 force fields for the physical relevant temperature range. Below 400 K these 
differences are within the errorbars, while for higher temperatures the ECEPP/3 energies 
are systematically lower by a small amount. Similar results were found for the specific 
heat C(T) as function of temperature shown in Fig. 7. Here we define 



Neither the position of its maximum nor its shape or heights vary much indicating that 
the transition between ordered and disordered states is of the same kind for both variants 
of the force field. Note that the position of the peak in specific heat corresponds to 
the increase of structures which belong to one of the main groups and the decrease of 
configurations which do not belong to any of the significant energy valleys. 

Conclusion 

We used the multicanonical approach to study variations of the energy landscape of Met- 
enkephalin under a change from the ECEPP/2 force field to ECEPP/3. In both cases, the 



C{(3) = 0- 



, 2 <E 2 >- 
5 



<E> 2 



(12) 
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low-energy local minima were sampled and classified into a small number of groups which 
correspond to essential valleys of the energy landscape. We studied their distribution and 
relative weight at various temperatures. While the energy landscape differs in detail and 
specially the global minimum is not the same, our physical results are little affected by 
changes between the two force fields. Our analysis demonstrates that at least for small 
peptides it is not sufficient to search only for the global minimum in potential energy since 
this conformation may be of limited significance at room temperature. Our work corrob- 
orates the well known fact that small flexible molecules exist at room temperature in an 
ensemble of low energy conformations (see for instance Ref. [|IJ and work cited therein) 
whose relative weight has to be determined. Finally, our results show that numerical sim- 
ulations of peptides are stable under small changes in the utilized force fields. Especially 
thermodynamic quantities are little changed. This is an important observation since all 
force fields rely on experimentally determined parameters and are therefore known only 
within certain errorbars. 
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TABLE CAPTIONS: 

TAB. 1: Point Charges for Terminal Groups of Met-Enkephalin. 

TAB. 2: Number of conformers and non-identical conformations collected after minimiza- 
tion for both choices of force fields. 

TAB. 3: Dihedral angles of our estimates of the global minima for ECEPP/2 and ECEPP/3. 
The peptide bond angles were fixed to oo = 180° in the simulations. 
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Table 1: 



Atom 


ECEPP/2 


ECEPP/3 


N 


-0.356 


-0.332 


H1,H2 


0.176 


0.076 


C 


0.450 


0.517 





-0.384 


-0.351 





-0.380 


-0.334 


H 


0.204 


0.235 
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Table 2: 





ECEPP/2 


ECEPP/3 


Minimized conformers E < —6 kcal/mol 


62142 


109863 


Non-identical minimized conformers 
(> 1 dihedral differs by > 1° ) 


1677 


1574 
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Table 3: 





ECEPP/2 


ECEPP/3 


1 Tyr 







-86.3 


-162.7 






T 


153.7 


-41.7 






UJ 


180.0 


180.0 






Yi 


-179.8 


-174.2 








-111.4 


-85.2 






AO 


145.3 


2.8 


2 Gly 




</> 

r 


-161.5 


65.8 






i> 


71.1 


-87.0 






UJ 


180.0 


180.0 


3 Gly 




<b 


64.1 


-157.3 






ib 

T 


-93.5 


34.9 






UJ 


180.0 


180.0 


4 Phe 




<b 


-81.7 


-158.8 






T 


-29.2 


159.5 






UJ 


180.0 


180.0 






Xi 


179.8 


52.4 






X2 


-100.0 


-96.0 


5 Met 







-80.7 


-82.4 








143.5 


134.1 






UJ 


180.0 


180.0 






Xi 


-65.1 


-66.1 






X2 


-179.2 


-179.6 






X3 


-179.3 


-179.9 






X4 


60.1 


60.1 


£ tot (ECEEP/2) [kca 


l/mol] 




-10.72 


-9.46 


^(ECEPP/3) [kca 


l/mol] 




-10.81 


-10.85 
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FIGURE CAPTIONS: 

FIG. 1: Timeseries of potential energy, obtained from a multicanonical simulation of 
200,000 MC sweeps, using a) the ECEPP/2 and b) the ECEPP/3 force field. 

FIG. 2: Timeseries of clusters, obtained from a multicanonical simulation of 200,000 MC 
sweeps, using a) the ECEPP/2 and b) the ECEPP/3 force field. 

FIG. 3: Most frequent low energy conformers for a) ECEPP/2 and b) ECEPP/3. Shown 
are also their Zimmermann codes, the frequency with which they appear in the multi- 
canonical simulation and the potential energy of their lowest energy exponent E min . 

FIG. 4: Groundstate configurations of a) ECEPP/2 and b) ECEPP/3. The plots were 
drawn using RASMOL. 

FIG. 5: Relative weight of the important groups of configurations in the case of a) 
ECEPP/2 and b) ECEPP/3. We have sum up all local minima which do not belong 
to one of the mayor groups in "Other". Likewise "No Minimum" contains all configu- 
rations which could not be minimized to a local minimum with potential energy below 
—6 kcal/mol (see text). 

FIG. 6: Average potential energy as a function of temperature for both choices of force- 
fields. 

FIG. 7: Average of specific heat as a function of temperature for both choices of force- 
fields. 
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